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It is shown that it is now possible to include van der Waals interactions via a nonempirical 
implementation of density functional theory to describe the correlation energy in electronic structure 
calculations on infinite systems of no particular symmetry. The vdW-DF functional [Phys. Rev. Lett. 
92, 246401 (2004)] is applied to the adsorption of benzene and naphthalene on an infinite sheet of 
graphite, as well as the binding between two graphite sheets. Comparison with recent thermal 
desorption data [Phys. Rev. B 69, 535406 (2004)] shows great promise for the vdW-DF method. 
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A recent study of the interaction of polycyclic aromatic 
hydrocarbon molecules (PAH's) with the basal plane of 
graphite Q] provides experimental benchmark data that 
constitute an ideal challenge for our recently proposed 
density functional (vdW-DF) which both includes van 
der Waals (vdW) interactions and all the benefits of the 
earlier state-of-the-art versions of the density-functional 
theory (DFT). Aiming at a better experimental charac- 
terization of the weak interlayer interactions in graphite, 
careful analysis of thermal-desorption kinetics yield acti- 
vation energies for benzene and PAH's at submonolayer 
coverages with explicit error bars |l( . Our calculated val- 
ues for the adsorption energy of benzene and naphtha- 
lene on graphene and for the weak interlayer interaction 
energy of graphene agree with the values deduced from 
experiment. From this we conclude that the vdW-DF 
is, indeed, very promising, and that it can be applied to 
systems that are neither periodic nor finite. This distin- 
guishes it from the various wave-function methods that 
are often applied to vdW complexes. 

Our method differs also from a newly published study 
of adenine on graphite [3j , which treats the vdW part of 
the correlation energy by a frequently used 0,|EIE0H)I3 
semi-empirical method. This method introduces empir- 
ical damping functions applied to an asymptotic attrac- 
tive 1/R 6 interaction assumed to occur between each 
pair of nuclei. At shorter distances this interaction is 
damped by a physically motivated, but arbitrary and 
varying functional form, which introduces one empirical 
parameter for every pair of atomic types in the complex. 
On the other hand our method (vdW-DF) for the cor- 
relation energy is completely free from empiricism, and 
although containing approximations, represents a first- 
principles density functional, which since the appearance 
of Ref . is applicable to arbitrary geometries, and which 
is seamless as two fragments merge into a single one. As 
discussed later, it has been applied to a number of phys- 
ical systems with promising results. The present appli- 
cation is particularly pertinent, however, as alternative 
first-principles methods for including vdW interactions 
are lacking for extended systems. 



Condensed matter is held together by several kinds of 
interatomic forces, including the ubiquitous vdW forces, 
which are particularly significant in sparse matter. For 
dense matter DFT has well proven its value, state-of-the- 
art versions of it giving values for ground-state properties 
of covalent molecules and hard materials close to exper- 
imental data. The key to success in DFT is the func- 
tional for exchange and correlation, and DFT calcula- 
tions today typically apply some flavor of the generalized- 
gradient approximation (GGA) 0, EH- For sparse 
matter, including soft matter, vdW complexes, poly- 
mers, biomolecules, and other abundant systems, how- 
ever, DFT in GGA performs badly. For instance, it give s 
unphysical results for the interlayer bond of graphite [12J > 
a canonical vdW case. The vdW interaction stems from 
truly nonlocal correlation and is contained in neither 
the local-density approximation (LDA) nor the semilocal 
GGA. We have subsequently proposed a density func- 
tional (vdW-DF) for general geometries Q that includes 
both vdW interactions and all GGA achievements. It 
has been applied to the benzene dimer 0, H^ ] and to 
complexes of benzene with phenol, toluene, fluoroben- 
zene, or benzonitrile 0]. All these previous applications 
have been to finite systems. We now apply it to single 
benzene and naphthalene molecules in interaction with 
an extended graphene sheet. This establishes the appli- 
cability of the vdW-DF functional to a more general class 
of systems. 

The derivation of the vdW-DF starts from the adia- 
batic connection formula and benefits from insights into 
the polarization properties of the inhomogeneous elec- 
tron gas. It involves several approximations, which sat- 
isfy known limits, sum rules, and invariances 0- To 
assess the validity of vdW-DF, internal criteria must 
be supplemented with external ones, that is, we must 
evaluate its performance in actual applications by com- 
parison with experimental findings. Systems where the 
vdW forces play a prominent role or dominate the in- 
teractions entirely are to be preferred. Being based on 
inhomogeneous-electron-gas polarization properties, the 
vdW-DF should work best on systems with similar po- 



2 



larization properties, that is, on delocalized systems with 
dense excitation spectra. Dimers of benzene and PAH's 
are such systems. Experimental data for these are scarce, 
and in practice the assessment has to be made by com- 
paring to accurate wave-function calculations 0, [H| , for 
which, however, the computational effort grows strongly 
with size, thus making benchmark results unavailable 
for really large systems. Further, to properly assess the 
vdW-DF method, extended systems should be addressed. 
The adsorption problem, with its extended substrate, is 
far beyond the possibilities of such wave-function meth- 
ods. The agreement between values calculated with 
vdW-DF and measured values of key quantities is demon- 
strated below. 

Recent thermal-desorption studies of the interaction of 
PAH's with the basal plane of graphite Q yield, through 
the analysis of desorption kinetics, activation energies 
of 0.50, 0.85, 1.40, and 2.1 eV for benzene, naphtha- 
lene, coronene, and ovalene, respectively, at submono- 
layer coverages. Using a force- field model the au- 
thors of Ref. ll| predict an exfoliation energy of 52 ±5 meV 
per/atom to peel a single graphene layer off a graphite 
surface. Another determination 16] involving collapsed 
nanotubes, which obtained 35j^Q meV/atom for a related 
quantity (the energy to separate two graphene sheets) , re- 
quired an equally tenuous theoretical model to extract a 
prediction. Here we apply vdW-DF to calculate the bind- 
ing energy of a benzene and a naphthalene molecule to a 
single graphene sheet, as well as the binding energy of two 
graphene sheets to each other. We also present evidence 
that second-layer interactions are small, which imply that 
our results also apply approximately to molecular des- 
orption or to graphene exfoliation from a basal graphite 
surface, as relevant for the experiment mentioned first 
above. 

The general geometry (gg) vdW-DF Q basically "cor- 
rects" the correlation part of the energy of a standard 
self-consistent (SC) DFT calculation, using the standard 
exchange from the GGA of the revPBE-flavor [l7|, cho- 
sen for its being closest to the Hartree-Fock exchange 
in some key applications |2j, [l3j. The vdW-DF energy 
functional is written 

-EvdW-DF = Egg & — -E-gga.c + -Blda.c + > (-0 

where the LDA correlation energy -Elda.c is substituted 
for the GGA correlation energy, E'qga.c, an d where the 
nonlocal correlation energy i?" 1 is added. The last three 
terms in are treated as a post-GGA perturbation, 
utilizing their low sensitivity to the choice of GGA SC 
electronic density. The nonlocal correlation energy is ex- 
pressed as 

Ef = \ J d 3 rd 3 r' ' n(r)<j>(r,r')n(r'\ (2) 

where <p(f,f') is a function depending on \f — r'\, the 
charge density n and its gradient at the positions r and 
r ', respectively, given in detail by Eq. (14) of Ref. 2. 



FIG. 1: The lateral configuration of the AB adsorption 
positions for the molecules benzene and naphthalene on a 
graphene sheet. Carbon atoms are shown in black, when 
belonging to the substrate and gray, when belonging to the 
adsorbate. Hydrogen atoms are shown in white. 



SC calculations are performed with the pseudo- 
potential-plane- wave-based DFT code DACAPO [lg, us- 
ing the revPBE functional to get GGA electron densities 
and total energies. The hexagonal cell used for the stan- 
dard part of the DFT calculations has the size (12.32 A, 
12.32 A, 26 A) for the benzene calculations and (14.783 
A, 12.32 A, 26 A) for the naphthalene calculations. We 
use a Monkhorst-Pack grid with (2, 2, 1) fc-points and a 
plane-wave cutoff energy of 450 eV. The in-plane struc- 
ture of the adsorbate is found through a monomer cal- 
culation, where the molecule is fully relaxed within a 
GGA-revPBE SC calculation. This results in an opti- 
mum structure with carbon-carbon and carbon- hydrogen 
bondlcngths in agreement with experimental data [19|. 
In the same way, the structure of the graphene sheet is 
obtained from an isolated calculation, resulting in an op- 
timum intralayer lattice constant of 2.46 A, in agreement 
with experiment |20j . Afterwards the in-plane structures 
are kept fixed while we vary the distance from the adsor- 
bate to the surface and map out an energy profile as a 
function of the separation. We place the molecule in an 
AB configuration (Fig.^l, which we expect to be energet- 
ically most favorable, as is the case for both benzene and 
naphthalene dimers, and for graphene layers. That is, we 
place the center of the benzene ring exactly above a car- 
bon atom in the graphene sheet, and for naphthalene, the 
molecule center-of-mass is positioned above the center-of- 
mass of the two graphite atoms below the naphthalene 
rings, as seen in Fig.n[2l]]. We know that for naphtha- 
lene dimers in the AB stacking, small shifts in the exact 
lateral position yield minor changes in the total energy, 
while the separation in the perpendicular direction is of 
much greater importance |22j| . and we believe the same 
should hold also for small lateral shifts with respect to 
the graphene layer in Fig. ^ 

For the adsorption position given above, full GGA cal- 
culations are performed, thus obtaining the first three 
terms in Eq. QJ. The last term, given by Eq. (J2J), re- 
quires some extra considerations, as (lateral) cell size 
becomes more important when evaluating the nonlocal 
vdW energy correction than for the standard calcula- 
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tions. For this purpose we enlarge the system size for the 
I?" 1 evaluation as follows. The SC electron density for the 
adsorbate-substrate system obtained with GGA-revPBE 
is used within the unit cell, while outside this unit cell the 
graphene electronic density is simulated by that of pure 
graphene, enabling the extension of the substrate to, in 
principle, arbitrary size |23| . Thus we can obtain the 
binding energy for benzene and naphthalene interacting 
with an increasingly larger circle of the graphene layer 
below the center of the adsorbed molecule. While 91% of 
the interaction is obtained for the original unit cell, 99% 
is obtained for a 10 A radius circle, and full convergence is 
reached for a 14 A circle. The same method is used for the 
calculation on the interactions of the two graphene layers. 
Thus calculation of E£ 1 allows an increase of system size 
without much increased computational cost compared to 
standard DFT calculations. 

Binding-energy results for two interacting graphene 
sheets in the AB stacking, which is the structure occur- 
ring in bulk graphite, are shown in Fig.0 The graphene- 
graphene binding-energy curve is somewhat deeper than 
that obtained with our earlier layered-functional version 
|12| . The general geometry version employs improvement 
made possible by an approximate expansion not made 
in the layered-functional version. Further discussion of 
the similarities and differences between the two methods 
is given separately [24(. The gg-vdW-DF curve deviates 
only slightly from the corresponding one in Ref. l25l where 
an averaging procedure was used. 

Our results for two graphene sheets in the AB stacking 
may be compared with the results of the experimental 
studies discussed above 0,^1- These are marked in Fig. 
|2]with diamonds at the layer separation in bulk graphite, 
which is 3.36 A .20} . We expect the separation of the AB 
graphene sheets to be similar, or only slightly larger. The 
energy error bars from the two experiments merge into 
one. 

While the exfoliation energy from a graphite basal 
plane is not the same as the energy to separate two 
graphene sheets, we may use the value of our curve at a 
second layer distance to estimate the change in our pre- 
diction if an exfoliation calculation were done instead. 
From the value of our curve of —2.5 meV at 7.2 A in 
Fig. [21 we can estimate that our method will give a 
48 meV/atom exfoliation energy, and a 50.5 meV/atom 
graphite cleavage energy, a 5.5% increase each [2(|. This 
conclusion requires the forces from subsequent layers to 
be negligible, which is certainly true for our functional, 
and probably also when the anomalous asymptotic effect 
of graphite's semimetallic nature j2]j is accounted for. 
For graphite, thermal vibrations also give a relevant con- 
tribution to the energy from the motion perpendicular to 
the sheets 7] . In any case the comparison of our curves 
with the points given by the two experimental groups is 
satisfying, despite the dependence of the latter on earlier 
theoretical methods that cannot be fully justified. 

Binding-energy curves are also calculated for benzene 
and naphthalene adsorbed on graphene, with configura- 
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FIG. 2: (Color online) Binding-energy curve for two graphene 
sheets in AB stacking. The calculated general geometry 
vdW-DF curve (solid line) indicates a binding energy of 45.5 
meV/atom at the equilibrium separation 3.6 A, which is close 
to the experimental energy estimates (diamonds with energy 
error bars; see text). The two vertical dashed lines mark the 
positions of the potential-energy minimum (3.6 A) and the ap- 
proximate point where a next-neighbor graphite layer would 
have been placed (7.2 A). 

tions as shown in Fig. We judge these calculations 
equivalent to those of adsorption on the basal plane of 
graphite, when a small correction for a second graphite 
layer (3.0% increase for benzene and 3.3% for naph- 
thalene) is taken into account, as discussed above for 
graphene-graphene. Figure |3| shows the calculated gg- 
vdW-DF potential-energy curves of adsorbed benzene 
and naphthalene, respectively, with their minima well 
within the ranges of measured binding-energy values 
0, l2q| . For benzene, we find a binding energy of 495 
meV (to be compared to the values 500 ± 80 meV 
and 480 meV J28J) at the equilibrium separation 3.6 A. 
The same separation is found optimal for naphthalene, 
with the adsorption energy 763 meV, which (particularly, 
when including the second-layer 3.3% contribution) can 
be compared to the experimental values 800 ±100 meV 
and 900 ± 70 meV resulting from slightly different analy- 
sis of the same experiment , summed up by the authors 
themselves by reporting the number 850 meV. The exper- 
imental values with error bars on the energy are shown as 
shaded regions in the figures. The separations were not 
measured in the experiment. However, it is reasonable 
to expect the separation to be similar to that in graphite 
also in these two cases, which is what we find. 

The analysis of the desorption experiments on ben- 
zene and naphthalene and the close agreement between 
experimental and theoretical results let us conclude that 
the vdW-DF is a very promising functional to account 
for the nonlocal vdW forces in these very typical cases. 
While the experimental evidence is not so direct in the 
graphene-graphene binding case, the values found are 
probably representative of what more direct experimen- 
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FIG. 3: (Color online) Benzene-graphene (left panel) and naphthalene-graphene (right panel) binding-energy curves. Like that 
for graphene-graphene, the general geometry vdW-DF calculations (solid lines) give results comparable to experiment (boxes); 
see text. The two vertical dashed lines mark the minima and the next-neighbor graphite layer, as in Figure |21 



tal methods would find, and give additional support to 
the vdW-DF method. 
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